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Abstract 

Regularization robust preconditioners for PDE-constrained opti¬ 
mization problems have been successfully developed. These methods, 
however, typically assume that observation data is available through¬ 
out the entire domain of the state equation. For many inverse prob¬ 
lems, this is an unrealistic assumption. In this paper we propose and 
analyze preconditioners for PDE-constrained optimization problems 
with limited observation data, e.g. observations are only available at 
the boundary of the solution domain. Our methods are robust with 
respect to both the regularization parameter and the mesh size. That 
is, the condition number of the preconditioned optimality system is 
uniformly bounded, independently of the size of these two parameters. 

We first consider a prototypical elliptic control problem and thereafter 
more general PDE-constrained optimization problems. Our theoretical 
findings are illuminated by several numerical results. 
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1 Introduction 


Consider the model problem: 


mm j^llu-dll^en)) + 

fll/ll!»(n)}’ 

(1) 

on a Lipschitz domain Cl C M n , subject to 



—An + u + / = 0 

in Cl, 

(2) 

du 


(3) 

Ihi = ° 

on dCl. 


This minimization task is similar to the standard example considered in 
PDE-constrained optimization. But instead of assuming that observation 
data is available everywhere in 0, we consider the case where observations 
are only given at the boundary dQ of 0, that is d E L 2 (d£i), see the first 
term in ([!]). For problems of the form ([!])- ([3]), in which 

-||u - d||^2( an ) +-||/||^2( n) (4) 

is replaced by 

-||n - cJ||^ 2 ( n ) +-||/|||2(Q) (5) 

very efficient preconditioners have been developed for the associated KKT 
system. In fact, by employing proper a-dependent scalings of the involved 
Hilbert spaces [TO], or by using a Schur complement approach [9], meth¬ 
ods that are robust with respect to the size of the regularization parame¬ 
ter a have been developed. More specifically, the condition number of the 
preconditioned optimality system is small and bounded independently of 
0<a<l and the mesh size h. This ensures good performance for suitable 
Krylov subspace methods, e.g. the minimum residual method (Minres), 
independently of both parameters. These techniques have been extended to 
handle time dependent problems [8j and PDE-constrained optimization with 
Stokes equations E2, but the rigorous analysis of a-independent bounds al¬ 
ways requires that observations are available throughout all of Cl. 

For cases with limited observations, for example with cost-functionals of 
the form 0. efficient preconditioners are also available for a rather large 
class of PDE-constrained optimization problems, see mm- But these tech¬ 
niques do not yield convergence rates, for the preconditioned KKT-system, 
that are completely robust with respect to the size of the regularization pa¬ 
rameter a. Instead, the number of preconditioned Minres iterations grows 
logarithmically with respect to the size of a -1 , as a —> 0: 

a + ftlogio ( a_1 ) • ( 6 ) 

1 In [6!, 7] it is proved that the number of needed preconditioned Minres iterations 


2 



According to the numerical experiments presented in [7], the size of b may 
become significant. More specifically, b G [5,50] for problems with simple 
elliptic state equations posed on rectangles. Thus, for small values of a, 
Minres may require rather many iterations to converge - even though the 
growth in iteration numbers is only logarithmically. 

In practice, observations are rarely available throughout the entire do¬ 
main of the state equation. On the contrary, the purpose of solving an 
inverse problem is typically to use data recorded at the surface of an ob¬ 
ject to compute internal properties of that object: Impedance tomography, 
the inverse problem of electrocardiography (ECG), computerized tomogra¬ 
phy (CT), etc. This fact, combined with the discussion above, motivate 
the need for further improving numerical methods for solving KKT systems 
arising in connection with PDE-constrained optimization. 

This paper is organized as follows. In the next section we derive the 
KKT system associated with the model problem 0-0- Our a robust pre¬ 
conditioner is presented in Section [3| along with a number of numerical 
experiments. Sections [4]{5] contain our analysis, and the method is general¬ 
ized in Sections [6](7j Section [8] provides a discussion of our findings, including 
their limitations. 

2 KKT system 

Consider the PDE 0 with the boundary condition ([3]). A solution u to 
this elliptic PDE, with source term / G P 2 (0), is known to have improved 
regularity, i.e. u G i7 1+s (0), for some s G [0,1], with s depending on 
the domain 0. In the remainder of this paper we assume that u has full 
regularity, i.e. u G H 2 (Q). This is known to hold if 0 is convex or if <90 is 
C 2 , see e.g. [2[ 3j. 

When solutions to ([2]) exhibit this improved regularity, we can write the 
problem on the non-standard variational form: Find u G H 2 ( O) such that 

{-Au + u,w) L 2 {n) + (f,w)u 2 (a) = o Vw G L 2 (0), (7) 


where 


H 2 { O) = <<f>€ H 2 ( 0) 


= 0 on <90 

an 


cannot grow faster than 

a + b [logic ( Q_1 )] ■ 

Furthermore, in [7[ it is explained why iterations counts of the kind © often will occur 
in practice. 
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equipped with the inner product 


(u,v) H 2 (n) = / V 2 ii : V 2 v + Vtt • Vu + uv dx 
■hi 


= / AuAv + Vu ■ Vu + uv dx. 

Jn 


( 8 ) 


Here V 2 rt denotes the Hessian of u, and the second identity is due to the 
boundary condition = 0 imposed on the space H 2 (Q). 

We will see below that, in order to design a regularization robust pre¬ 
conditioner for ([I])-([ 3 ]), it is convenient to express the state equation in the 
form ([7]), instead of employing integration by parts/Green’s formula to write 
it on the standard self-adjoint form. 


2.1 Optimality system 

We may express Q-Q in the form: 

41 2 <“> + W 

subject to 

(-Au + u,w) L 2 ^ + (f,w) L 2 (n) = 0 VmGt 2 (0). (10) 


The associated Lagrangian reads 


1, 


a , 


C(f,u,w) = ~\\u - d\\l 2(gn) + -||/|| L 2 (n) + (/ - A u + u,w) L 2 {a) , 

with / E L 2 (H), u E H 2 (£l) and w E L 2 (H). From the first order optimality 
conditions 

^ = 0 ^ = 0 ^=0 
df ’ du ’ dw 

we obtain the optimality system: Determine ( f,u,w ) E L 2 (Q) x H 2 ( H) x 
L 2 (Q) such that 


(u 


a(/,^)L2(f 2 ) + 

d, ( t ) ')L 2 {dQ.) + (—A4> + (f>, w) L 2^ 

{h€)L 2 (Q) + ( —Ari + U, £)l 2 (T2) 


= 0 e L 2 (n), 

(11) 

= 0 V</> E H 2 (n), 

(12) 

= 0 V£ E L 2 (Q). 

(13) 


3 Numerical experiments 

Prior to analyzing our model problem, we will consider some numerical 
experiments. Discretization of yields an algebraic system of the 
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form 


aM 

0 

M " 


/ 


0 

0 

Mq 



u 

= 

M d d 

M 

A 

0 


w 


0 


( 14 ) 


A a 


where 


• M is a mass matrix, 

• Mq is a mass matrix associated with the boundary dQ of Q. 

• A is a matrix that arise upon discretization of the operator (1 — A). 
Since we write the state equation on a non self-adjoint form, A will not 
be the usual sum of the stiffness and mass matrices. Instead, equation 
([7]) is discretized with subspaces of H 2 ( fl) and L 2 (Vt). 

In the current numerical experiments, we employ the Bogner-Fox-Schmit 
(BFS) rectangle for discretizing the state variable u E H 2 (Q). That is, the 
finite element field consists of bicubic polynomials that are continuous, have 
continuous first order derivatives and mixed second order derivatives at each 
vertex of the mesh. BFS elements are C 1 on rectangles and therefore H 2 - 
conforming. The control / and Lagrange multiplier w are discretized with 
discontinuous bicubic elements. 


We propose to precondition (14) with the block-diagonal matrix 


B n = 


aM 0 0 

0 aR + Mg 0 
0 0 


-1 -1 


-M 


(15) 


where R results from a discretization of the bilinear form £>(•, •) on H 2 (Q): 

b(u,v) = (u,v)h 2 (q) + / Vu-Vvdx. (16) 

Jn 

In the experiments presented below, we used this bilinear form to construct 
a multigrid approximation of (aR + Mq)^ 1 . 


Remark 


The bilinear form (16) is equivalent to the inner product on iL 2 (fl). The 


additional term stems from our choice of implementing a multigrid algorithm 
for the bilinear form associated with the operator (A — l) 2 = A 2 — 2A + 1. 
Indeed, the bilinear form ab( •, •) + (■, • )l 2 (QQ) can be seen to coincide with 
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the variational form associated with the fourth order problem 


a(A — l) 2 u 
du 

dn 

dAu 
a- 


<9n 


= / 

in II, 

= 0 

on d£l, 

= u 

on dQ. 


To limit the technical complexity of the implementation, we considered the 
problem ([Tj)- ([3]) on the unit square in two dimensions. The experiments were 
implemented in Python and SciPy. The meshes were uniform rectangular, 
with the coarsest level for the multigrid solver consisting of 8 x 8 rectangles. 
Figure [I] shows an example of a solution of the optimality system (14). 


3.1 Eigenvalues 


Let us first consider the exact preconditioner B a defined in (15). If B a is a 


good preconditioner for the discrete optimality system (14), then the spectral 
condition number of B a A a should be small and bounded, independently 
of the size of both the regularization parameter a and the discretization 
parameter h. 

The eigenvalues of this preconditioned system were computed by solving 
the generalized eigenvalue problem 


A a x = A B^x. 

We found that the absolute value of the eigenvalues A were bounded, with 

0.445 < |A| < 1.809, 

uniformly in a G {1,10 -1 ,..., 10 10 } and h E {2“ 2 ,..., 2 -5 }. This yields 
a uniform condition number k(B a A a ) ~ 4.05. The spectra of the precondi¬ 
tioned systems are pictured in Figure [2] for some choices of a. The spectra 
are clearly divided into three bounded intervals, and the eigenvalues are 
more clustered for a ~ 1 and for very small a. 

3.2 Efficient preconditioning 

In practice, the action of B a is replaced with a less computationally ex¬ 
pensive operation B a . Note that B a has a block structure, and that com¬ 
putationally efficient approximations can be constructed for the individual 
blocks. Specifically, B a is constructed by employing 
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Observation data Optimal state 



(a) Observation data d. The forward model (b) Computed optimal state u based on the 

was solved for the control shown in |d]), but observation data shown in ja|) 

only the boundary values can be observed. 


Optimal control Generating control 




(c) Computed optimal control / based on the 
observation data in ja|) 


(d) The “true” control function used to gen¬ 
erate the observation data in 13 - 


Figure 1: An example of a solution of (14). The observation data d was 


generated with the forward model, using the “true” control 4x(l — x) + y 
shown in panel ([d]). Solutions to the unregularized problem are non-unique, 
and the generating control cannot be (exactly) recovered. The figures were 
generated with mesh parameter h = 1/128 and regularization parameter 
a = 10 -6 . 
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0.445 

0.000 


Sorted eigenvalues for a=10° 


4000 6000 

k 



(a) a = 1 (b) a = 10 4 




(c) a = 10“ 6 (d) a = 10“ 10 

Figure 2: Spectrum of B a A a for different regularization parameters a. The 
discretization parameter was h = 2 -4 for all figures. 
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Iterations 

1 

2 

3 

(h = 2"«) 

1.931 

1.303 

1.126 


Table 1: Condition numbers of 
M preconditioned with symmetric 
Gauss-Seidel iterations. 


a\h 

2 -4 

2 -6 

2 -s 

1 

1.130 

1.136 

1.140 

10“ 4 

1.129 

1.135 

1.139 

10 -8 

1.237 

1.150 

1.149 

10~ 12 

1.252 

1.259 

1.253 


Table 2: Estimated condition 
numbers of aR + Mg precondi¬ 
tioned with one V-cycle multigrid 
iteration. 


• 1 multigrid V-cycle for the (2,2) block of B a , containing a symmetric 
4x4 block Gauss-Seidel smoother where the blocks contain the matrix 
entries corresponding to all degrees of freedom associate with a vertex 
in the mesh (see EH for a theoretical analysis of the method). 

• 2 symmetric Gauss-Seidel iterations for the (1,1) and (3,3) blocks. 

We estimated condition numbers of the individual blocks of R” 1 precon¬ 
ditioned with their respective approximations. The results are reported in 
Tables [l] and [2j A slight deterioration in the performance of the multigrid 
cycle can be seen for very small values of a > 0. 

3.3 Iteration numbers 

To verify that also B a is an effective preconditioner for A a , we applied the 
Minres scheme to the system 

B a A a x = B a b. 

For the results presented in Table [3j the Minres iteration process was 
stopped as soon as 

(rk,B a r k ) _ (■ AgX k - b, B a {A a x k - b }) < ^ ^ 

(r 0 , B a r 0 ) {A a x o - b, B a {A a x 0 - b }) 

which is the standard termination criterion for the preconditioned MlNRES 
scheme, provided that the preconditioner is SPD. A random initial guess xo 
was used, and the tolerance was set to e = ICG 12 . 

4 Analysis of the KKT system 

Recall that our optimality system reads: 

a(f, VOlRo) + Wb w) L 2 {n) =0 € L 2 (n), 

( u — d, 4>)L 2 (an) + + 0; w )l 2 ( n) =0 W> £ H 2 (Cl), 

(/) 0l 2 (O) + (—Au + u, 0i 2 (a) =0 e L~(Q), 
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a\h 

2 -4 

2" 5 

2 -6 

2 -7 

1 

53(4.33) 

53(4.36) 

53(4.36) 

53(4.36) 

10” 1 

57(4.31) 

57(4.34) 

57(4.35) 

57(4.35) 

10 -2 

75(4.31) 

72(4.34) 

70(4.35) 

68(4.35) 

10" 3 

79(4.31) 

79(4.34) 

77(4.35) 

73(4.35) 

10 -4 

81(4.30) 

81(4.33) 

79(4.35) 

77(4.35) 

10 -5 

82(4.33) 

81(4.33) 

79(4.35) 

79(4.35) 

10" 6 

81(4.35) 

79(4.36) 

79(4.35) 

81(4.35) 

10- 7 

70(4.35) 

81(4.37) 

81(4.36) 

79(4.35) 

10“ 8 

62(4.36) 

70(4.36) 

79(4.36) 

81(4.36) 

10“ 9 

62(4.36) 

64(4.37) 

68(4.37) 

78(4.36) 

10- 10 

62(4.36) 

63(4.36) 

64(4.37) 

67(4.37) 


Table 3: Number of preconditioned Minres iterations needed to solve the 
optimality system to a relative error tolerance e = 10 -12 . Estimated condi¬ 
tion numbers in parentheses, computed from conjugate gradient iterations 
on the normal equations for the preconditioned optimality system. 


with unknowns / E L 2 {Q), u E 77 2 (fl) and w E L 2 (0). We may write this 
KKT system in the form: 

Determine ( f,u,w ) E L 2 (fl) x H 2 (£l) x L 2 (D) such that 


1 

P 

o 

_ 1 


' / ' 


0 

0 M d A 1 


u 

= 

M d d 

L M A 0 


w 


0 


-Ac* 


where 


(18) 


M : 

: L 2 (fl) - 

->• L 2 m\ 

f ^ (/> • )l 2 {CL)i 

(19) 

M d : 

: H 2 (Q) 

->• H 2 (fi)', 

(u, ■ )l2 ( 9 q), 

(20) 

Mq : 

b-i 

to 

5 

ff 2 (n)', 

d | -t (d, ■ )L 2 {dCl): 

(21) 

A : 

: H 2 (Q) 

L 2 (Q)', 

u (-A u + u, -)l 2 (Q) 5 

(22) 


and the notation ” / ” is used to denote dual operators and dual spaces. In the 
rest of this paper, the symbols M, Mg and A will represent the mappings 

), respectively, and not (the associated) matrices, 
(We believe that this mild ambiguity improves 


defined in (fl9|) , pOl) and 


as was the case in Section D 
the readability of the present text). 

By using standard techniques for saddle point problems, one can show 
that the system (18) satisfies the Brezzi conditions [I], provided that a > 0. 
Therefore, for every a > 0, this set of equations has a unique solution. 
Nevertheless, if the standard norms of L 2 (D) and H 2 (Q) are employed in the 
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analysis, then the constants in the Brezzi conditions will depend on a. More 
specifically, the constant in the coercivity condition will be of order O(a), 
and thus becomes very small for 0 < a < 1. This property is consistent 
with the ill posed nature of Q-Q for a = 0, and makes it difficult to design 
a robust preconditioners for the algebraic system associated with (18). 


Similar to the approach used in [51 El HD] , we will now introduce weighted 
Hilbert spaces. The weights are constructed such that the constants appear¬ 
ing in the Brezzi conditions are independent of a. Thereafter, in Section [5] 
we will show how these scaled Hilbert spaces can be combined with simple 
maps to design a robust preconditioners for our model problem. 


4.1 Weighted norms 

Consider the a-weighted norms: 


/ L 2 (n) — a \\f li 2 (n)> 

(23) 

\\ u H 2 (n) = a W u // 2 (0) + M £2(90), 

(24) 

; ll^_i(n) = “IMI l 2 (Q)> 

(25) 


applied to the control /, the state u and the dual/Lagrange-multiplier w, 
respectively. Note that these norms become “meaningless” for a = 0, but 
are well defined for positive a. 


4.2 Brezzi conditions 

We will now analyze the properties of 

A a : L 2 (Q) x H 2 (Q) x L 2 a ^(Q) x ^(H)' x L^H)', 


defined in (18). More specifically, we will show that the Brezzi conditions are 
satisfied with constants that do not depend on the size of the regularization 
parameter a > 0. Note that we use the scaled Hilbert norms (23)-(25). 

Lemma 4.1 For all a > 0, the following “inf-sup” condition holds: 

(/, w )l 2 {0) + (-Au + u, w) L 2 (n) 


inf 


eL 2 _ 1 (n)(/,u)ei|(Q)xif2(n) \\(f^ u )\\Ll(a)xH^)W w \\L 2 ^(fi) 


> 1. 


Proof 

Note that L^(Q) and F 2 a _ 1 (yt) contain the same functions, provided that 
a > 0. Let w G L 2 a _ l (12) be arbitrary. By choosing f = w and u = 0 we 
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find that 


sup 


(/» w)l2(Q) + (-■ Au + ™)l 2 (Q) 


> 


(w,w) L 2 (n) 


(/,«)e-L2(n) x /t 2 (Q) IK/ ! ' u )llL2 ( n) xH 2(Q)||u;|| L 2_ i(n) IK^)0)llL 2 (n)xi7 2 (n)ll' u; llL2_ 1 (a) 




IIl 2 (Q) 


V^Mli 2 ^)^) MMIl 2 ^) 

= 1. 


Since u> E -L 2 „i(ri) was arbitrary, this completes the proof. 


Expressed in terms of the operators that constitute A a , Lemma 4.1 takes 
the form 


inf 


sup 


(Mf, w) + ( Au , w) 


sl 2 _ 1 (O) (/,«) ei 2 (a) x h 2 (Q) 11 (/,«) 11 lz (n) x H 2 (n) 11 w 11 L 2 _ i (Q) 


> 1, 


see (19) and (22). 

Recall that we decided to write our state equation ([2|-([3]) on the non¬ 
standard variational form ([7]) . Throughout this paper we assume that prob¬ 
lem <©-©> admits a unique solution u E H 2 (Q) for every / E L 2 (fl), and 
that 

IMIff 2 (n) < c iII/IIl 2 (q)- (26) 

This assumption is valid if fl is convex or if fl has a C 2 boundary, see e.g. 
Inequality ( |26[ ) is a key ingredient of the proof of our next lemma. 

Lemma 4.2 There exists a constant C 2 , which is independent of a > 0, 
such that 


a ll/llL 2 (n) + IMIl 2 (<90) ^ c 2 {^W/Wl 2 ^) + a \\ u \\H 2 (CL) + ll u 

= c 21| (/, u )\\l 2 (Q)xH 2 (CI) 
for all ( f,u ) E L 2 ( Q) x H 2 (Q) such that 

(/> 4>)L 2 (n) + (~ Au + u, 4 >)l 2 {q.) =6 V0 E L 2 (p). 


(27) 


Proof 


If (f,u) satisfies (27), then 

IMIir 2 (n) - c iII/IIl 2 (q)) 

see the discussion of (26). Let 0 = (1 + c 2 ) -1 e(O, 1), and it follows that 

“ll/lli 2 ^) + IMIl 2 ((9 q) > CK0 ||/||^2(Q) + a —^r\\ u \?H 2 {Q) + \\ u \\l 2 (dn) 


> 


1 


1 + cf 


2 ( a ll/ll22(Q) + a \\ u \\ 2 H 2 (Q.) + ll^" 2 
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This result may also be written in the form 


aM 0 1 1/1 l"/l\ n m 2 

0 Mg\ ’ [«J/ - ^(n)xfl2(n) 

for all (/, u) € L^(Q) x H 2 (Q) satisfying 

Mf + Au = 0, 


where M, Mg and A are the operators defined in (19), (20) and (22), re¬ 
spectively. 


4.3 Boundedness 

Having established that the Brezzi conditions hold, with constants that are 
independent of a , we next explore the boundedness of A a . 


Lemma 4.3 


/ 

aM 

0 

f 


-if 

\ 

\ 

0 

Mg 

u 

5 


/ 



for all (/,«), (?M) G L 2 a (n) x 


Proof 


Recall the definitions (19) and (20) of M and Mg, respectively. Since 

IKMII L2(o)x//2(Q) - \J “II/IIlro) + a IMIjf 2 (n) + IMli 2 (dn)> 

we find, by employing the Cauchy-Schwarz inequality, that 


/ 

aM 

0 

f 


-if 

\ 

0 

Mg 

u 

1 

<\> 


— |a(/,^)L2(n) + ( u > 4>)L 2 (an)\ 

^ 11/11x2(0)11^11x2(0) + IMIx 2 (,9Q) II 11x2(90) 

^ ^sJWfWli (0)11^11x2(0) + IMIx2(do)ll/ , llx2(dn) 

< V / 2||(/,u)|| L 2 (Q) X Hi (O) 11 11X 2 (0) x 2 (O) • 


Lemma 4.4 


[M H] 




< V ^ IK /, u )|| L 2 ( 0 ) xJX 2 ( 0 ) 11 ^ 11 x 2 _ i ( 0 ) 


for all (f,u) € L£(fi) x #2 (n)> w g L 2 (n)i 
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Proof 

Again, we note that 


II (/; ^) llz,2 (f2)xj/2(Q) — \J a \\f\\L 2 {u) + a IMI_fl- 2 (n) + IMIz, 2 (dn)> 


l”'ll^-i(n ) 


\w 


IIl 2 (0)- 


From the definitions of M and A, see (19) and (22), and the Cauchy-Schwarz 
inequality, it follows that 


[M A] 


, w 


= | (Mf,w) + (Au,w) | 

= | (/, w)l*(si) + (-■Au + U, w) L 2 (n) I 

- (||/IIl 2 (o) + II Au IIl 2 (o) + II^IIl 2 (o)) 

— v/ 3|l(/!' u )llL 2 (n)x_ff 2 (n)ll u, llL 2 _ 1 (n)- 


ufi 7-2 


i-lM 


For the last equality, recall from © that j|Au|| L2(n) = ||V 2 u|| L2(n) < 

IMI.ff 2 (n) f° r a ll u £ H 2 (Q). 


4.4 Isomorphism 

We have verified that the Brezzi conditions hold, and that A a is a bounded 
operator. Moreover, all constants appearing in the inequalities expressing 
these properties are independent of the regularization parameter a > 0. Let 

v = L 2 a (n) x H 2 a (n) x ( 28 ) 

V = L 2 a (nyxH 2 a (nyxL 2 a ^(ny. (29) 


Theorem 4.5 The operator A a , defined in (18), is bounded and continu¬ 
ously invertible for a > 0 in the sense that for all nonzero x E V, 

(A q x, y ) 


c < sup 

o^2/ev 


<C, 


(30) 


/or some positive constants c and C that are independent of a > 0. In 
particular, 

II ll£(V',V) — C aTl d 11 11 J C(V,V") — C- 


Proof 


This result follows from Lemma 4.1, Lemma 14.21 Lemma 4.3, Lemma 4.4 


and Brezzi theory for saddle point problems, see [Tj. 
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4.5 Estimates for the discretized problem 


The stability properties (30) is not necessarily inherited by discretizations. 
However, the structure used to prove the so-called “inf-sup condition” in 


Lemma |4.1| is preserved in the discrete system provided that the same dis¬ 
cretization is employed for the control and the Lagrange multiplier. Further¬ 
more, the boundedness properties, Lemma 4.3 and Lemma 4.4, certainly also 
hold for conforming discretizations. 

It remains to adress the coercivity condition, Lemma 4.2, for the dis¬ 
cretized problem. We consider finite dimensional subspaces Uh. C U = 
H 2 {Q.) and Wh C W = L' 2 (Q). For certain choices of Uh and Wh, the 
estimate of Lemma 4.2 carries over to the finite-dimensional setting. 


Lemma 4.6 Assume Uh C U and Wh C W , such that (1 — A )Uh C Wh- 
Then 

a 11 All 1,2(0) + \\ u h\\L2(dil) - c 2||(//ii^/i)lli 2 (Q)xH 2 (f2) (31) 

for all ( fh,Uh ) G I Vh x Uh such that 


(, fh, <A)l 2 (0) + ( u h - A U h , 4>h) L 2 (Q) = 0 ^ <fh G W h - (32) 


Proof 


Assume that (1 — A )Uh C Wh, and that (32) holds for ( fh, u h ) G Wh x Uh- 
Then fh + (1 — A )uh G Wh, and (32) implies fh + (1 — A )uh = 0. Therefore, 
( fh, u h .) satisfies (|27|) and the estimate (31) follows from Lemma 4.2 


If the discretization is chosen such that Lemma 4.6 is satisfied, then the 


estimates (|30|) carries over to discretized system. More precisely, we have 

and 


\\-A-a.,h\\c(y h ,V' h ) — ll“^«ll£(V,V')’ 


IIA*,J>C(VLV h ) - 11A 


-il 
a \ 


£(V',V)> 


(33) 

where Vh = Wh X Uh x Wh C V, equipped with the inner prdocut of V, 
and A at h is discrete counterpart to A a , defined by setting {A a .hXh,Uh) = 
(■ A a x h ,yh ) for all x h ,Vh G V h - 

If the state is discretized with C' 1 -conforming bicubic Bogner-Fox-Schmit 
rectangles, as in Section [3j then Lemma 4.6 is satisfied if the control and 
Lagrange multiplier is discretized with discontinuous bicubic elements on 
the same mesh. For triangular meshes, one could choose Argyris triangles 
for the state variable and piecewise quintic polynomials for the control and 
Lagrange multiplier variables. 


We remark that Lemma 4.6 provides a sufficient, but not necessary cri¬ 
terion for stability of the discrete problem, and usually may imply far more 
degrees of freedom in the discrete space Wh C W than is actually needed. 
The usefulness of Lemma 4.6 is that the estimates (33) can, in principle, 


always be obtained by choosing a sufficiently large space for the control and 
Lagrange multiplier. 
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5 Preconditioning 

The linear problem ( |18[ ) is of the form 

Ax = b. (34) 


where x is sought in a Hilbert space V, the right hand side b is in the dual 
space V', and A is a self-adjoint continuous mapping of V onto V. Iterative 
methods for linear problems are most often formulated for operators map¬ 
ping V into itself, and can not be directly applied to the linear system (34), 
as described in [5]. If we want to apply such methods to (34), then we need 
to introduce a continuous operator mapping V isomorphically back onto V. 
More precisely, if we have a continuous operator 

B : V ->• V, 


then M. = BA : V —> V is continuous and has the desired mapping prop¬ 
erties, and if B is an isomorphism, the solutions to (34) coincides with the 
solutions to the problem 


Mx = BAx = Bb. (35) 

In this paper we shall consider B G £(V,V 7 ) a preconditioner if B is 
self-adjoint and positive definite. This implies that B is self-adjoint and 
positive definite as well, and hence B defines an inner product on V by 
setting 

(x, y) = (B~ x x, y), x,yeV. (36) 

This inner product has the crucial property of making M self-adjoint, in 
the sense that 


(Mx, y) = (Ax, y) = (Ay,x) = (My,x). (37) 


Conversely, given any inner product on (•, •) on V, the Riesz-Frechet 
theorem provides a self-adjoint positive definite isomorphism B : V —> V 
such that ( p6| ) and (37) hold, and we say that B is the Riesz operator 
induced by (•, •). This establishes a one-to-one correspondence between 
preconditioners and Riesz operators on V 7 . Since the Riesz operator is an 
isometric isomorphism, the operator norm of BA coincides with the opera¬ 
tor norm of A. We formulate this well-known fact here in a lemma for the 
sake of self-containedness. We refer to 1313 for a more in-depth discussion 
of preconditioning and its relation to Riesz operators. 


Lemma 5.1 Let V be a Hilbert space, and let A : V —>• V be a self-adjoint 
isomorphism, and assume that B is the Riesz operator induced by the inner 
product on V, or equivalently, that the inner product on V is defined by the 
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self-adjoint positive definite isomorphism B 1 : V —> V'. Then BA : V —> V 
is an isomorphism, self-adjoint in the inner product on V, with 

ll^llz:(v,v) = ll-4|l.c(v,v') an d ll(^*A) 1 llz:(v,v) = ll-^ 1 ll£(v , ,v)- 

In particular, the condition number of BA is given by 

k(BA) = ||-4 1 ||£(v , ,v)ll'4|l£(v,v / )' 


Proof 


Since A is self-adjoint, M = BA is self-adjoint with respect to the inner 
product on V. From the Riesz-Frechet theorem we have ||-4x|| v , = ||£14x|| = 
||.A/fx||, and we obtain following identity for the operator norm of Ai. 


\\Mx\\ v ||Ac|| v , 

£(V,V) = SU P II. .11 = SU P II..II 


= sup sup 


{-Ax, v) 


c{vy)- 


x^O y ^=0 |FHvll2/llv 

A similar identity is obtained for the norm of the inverse operator, 

-in _ll-^^llv 


II M 


\C{V,V) = SU P 

x^O 





-i 


= inf sup 


{Ax, y) 


l x^O IFlIvllJ/llv 


= M 


-ii 


l£(V',V)' 


We say that a preconditioner B a for A a is robust with respect to the param¬ 
eter a if Hi{B a A a ) is bounded uniformly in a. The significance of Lemma |5.1| 
is that such a robust preconditioner can be found by identifying (parameter- 
dependent) norms in which A a and Ajj 1 are both uniformly bounded. 


5.1 Parameter-robust minimum residual method 


In Section [4] stability of A a was shown in the a-dependent norms defined in 


The preconditioner provided by Lemma 5.1 is the Riesz operator 
induced by the weighted norms. This operator B a ■ V —> V takes the form 


B 


a 


aM 0 0 

0 aR + Mq 0 
0 0 -M 

a 


(38) 
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where R : H 2 {Lt) —> H 2 (Ll)' is the operator induced by the H 2 (Q) inner 
product, i.e. (Ru,v) = (u,v) H 2 ^y 

Since A a is self-adjoint, the preconditioned operator B a A a ■ V —> V 
is self-adjoint in the inner product on V. Consequently we can apply the 
minimum residual method (MlNRES) to the problem 

B a A a x = Bab. 


Theorem 5.2 Let A a be the operator defined in (18) and B a the operator 
defined in (38). Then there exists an upper bound, independent of a, for the 
convergence rate of Minres applied to the preconditioned system 


B„A a x = B„b. 


In particular there exists an upper bound, independent of a, for the number 
of iterations needed to reach the stopping criterion ©• 


Proof 

A crude upper bound for the convergence rate (more precisely, the two-step 
convergence rate) of Minres is given by 


II B a Aa{x 


X2m) ||y < 


1 — K 
1 + K 


\\B a A a (x 


*o)|lv 


where n = K,{B a A a ) is the condition number of B a A a , see e.g. [5j. From 
Lemma |5.1| and pO] ) we determine that k is bounded independently of a , 
with 

K= | \(B a A a ) 1 ||£(v,v)II^«“^“II£(v,v) 

= 1 ll/:(v , ,v) ll^-« ll/ztv.vo (39) 

< c- x C. 


In practical applications, the operator B a will be replaced with a less compu¬ 
tationally expensive approximation B a ■ Ideally B a will be spectrally equiv¬ 
alent to B a , in the sense that the condition number of B a B~ l is bounded, 
independently of a. Then the preconditioned system reads 

B a A a x = Bab , 

and the upper bound for the convergence rate is determined by the condi¬ 
tioned number K(B a A Q ) < K{B a B- l )K{B a A~ l ). 


18 





Remark 


In this paper we only consider the minimum residual method, and we there¬ 
fore require that the preconditioner is self-adjoint and positive definite. More 
generally, if other Krylov subspace methods are to be applied to (18), then 


preconditioners lacking symmetry or definiteness may be considered. 

We mention in particular that a preconditioned conjugate gradient method 
for problems similar to (18) was proposed in [10] . based on a clever choice 
of inner product. 


6 Generalization 

Is our technique applicable to other problems than Q-Q? We will now 
briefly explore this issue, and show that the preconditioning scheme derived 
above yields a robust methods for a class of problems. 

The scaling (23)-(25) was also investigated in [6], but for a family of 
abstract problems posed in terms of Hilbert spaces. More specifically, for 
general PDE-constrained optimization problems, subject to Tikhonov regu¬ 
larization, and with linear state equations. But in |6j no assumptions about 
the control, state or observation spaces were made, except that they were 
Hilbert spaces. Under these circumstances, it was proved that the coerciv- 
ity and the boundedness, of the operator associated with the KKT system, 
hold with a-independent constants. Nevertheless, in this general setting, 
the inf-sup condition involved an a-dependent constant, which, eventually, 
yielded theoretical iteration bounds of order 0([log (a -1 )] 2 ) for Minres. 

In the present paper we were able to prove an a-robust inf-sup condition 
for the model problem ©-©• This is possible because both the control / 
and the dual/Lagrange-multiplier w belong to L 2 (H). From a more general 
perspective, it turns out that this is the property that must be fulfilled in 
order for our approach to be successful: The control space and the dual 
space, associated with the state equation, must coincide. This will usually 
lead to additional regularity requirements for the state space. 

Motivated by this discussion, let us consider an abstract problem of the 
form: 

(40) 


mm 

few,ueu 


i^WTu-d\\ 2 0 + ^a 


subject to 


{Au, w) + (/, w)w = 0, Vu; E W. 


(41) 


Here, 


W is the dual and control space, 
U is the state space, 


O is the observation space, 
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• W, U and O are Hilbert spaces. 

Let us assume that 

(Al) A : U —> W' is a continuous linear operator with closed range. In 
particular, there is a constant c\ such that for all u E U, 

\\ u \\u/KevA = inf 11“ - u\\u ^ c A\ Au \\w 
itGKer A 

(A2) T :[/—>• O is linear and bounded, and invertible on the kernel of A. 
That is, there is a constant C 2 such that for all u E Ker A, 

Nlcr < c 2\\Tu\\ 0 . 


It then follows that the KKT system associated with (140j)-( 41) is well-posed 


for every a > 0: Determine (/, u, w) E W x U x W such that 


aM 

0 

M r 


r/i 


' 0 ' 

0 

K 

A 1 


u 

= 

kd 

M 

A 

0 


w 


_ o _ 


(42) 


=4 q 


where 


M :W - 


f ^ if, ■ )w, 

(43) 

K :U - 


u (Tu, T ■ ) 0 , 

(44) 

>7 

O 

1 

>u’, 

d e->• (d, T ■ )o, 

(45) 


Note that, compared with (14), the boundary observation matrix Mg has 


been replaced with the general observation operator K in (42). 
We introduce scaled norms as follows. 


2 

w a 

2 


2 

W’ 
2 


= a | 

u \\u a = a ll^- u lliy' + Iloilo, 

2 1 n 1.2 


|W|Im/ -- 1 “ a'' W " w ' 


We first show that 


| jj is indeed a norm on U when assumptions 
and (A2) hold. It suffices to show that || • |I{/ Q is a norm equivalent to 


(Al) 


\u 


when a = 1. We have 

Iloilo + ll^llw' — (ll^ll C(U,0) + ll-^ll C(U,W')) II “11(7’ 
and letting n denote the orthogonal projection of U onto Ker A , 


u 


(46) 


M|[/ < W^uWjj + 11 It — 7TU 11 jj 


< c 2 ||T7ru|| 0 + ||u - Trull y 

< c 2 ||Tu|| 0 + (l + c 2 ||T|| £{C7)0) )||u-7 r u||j 7 

< c 2 ||Tu|| 0 + ci(l + C 2 H^H,c(c/,o)) ll-^llw' - 


( 47 ) 
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Here the last inequality follows from \\u — 7ru||^ = inf SgKer a\\ u ~ ^Wu aR d 
assumption |(A1)| 

We set V = W a xU a x W a -i . As in Section[4j A a : V —> V can be shown 
to be an isomorphism, with parameter-independent estimates obtained in 
the weighted norms. 

Theorem 6.1 There exists positive constants c and C, independent of a, 
such that for all nonzero x 6 V, 


( A a x , y) 

C < SUP 77-Jj-Jj- n - 

oI h 


< C. 


(48) 


ivi 


V 


We omit the full proof, which is analogous to that of Theorem 4.5 The 
crucial part is the “inf-sup condition” of Lemma 4.1. which is easily shown 
to hold in the abstract setting: 


(. f,u)ew a xu a W(fi u )\\w a xU a \\' w \\w n -1 ll( u; > 0 )llw a x! 7 a ll' u; llm _1 


(/> uj) w + {Au,w} 


> 


(w,w] 


w 


= 1. 


The coercivity condition of Lemma 4.2 naturally holds in the prescribed 
norm on U a , since for ( f,u)GWxU such that Au = Mf , 


<*tt.T\\w + Ml = wll/llvy + ^WAufw + \\Tu\\ 2 0 > - 


2 i II || 2 

w a + IMIt / 0 


Note that the weighted norm now depends on A, and as consequence, 
the estimates become A-independent. In fact, we obtain bounds for the 
constants c and C which are independent of a as well as the operators 
appearing in (40)-(41). This is postponed to the next section, where sharp 
estimates are obtained for (48). 

With the estimates (48), Lemma 5.1 provides a preconditioner for the 
operator A a , given as 


B n = 


aM 0 

0 aA’M~ 1 A + K 
0 0 


0 
0 

-M 

a J 


-i -1 


(49) 


The condition number of B a A a will be bounded independently of a. It is, 
however, not clear how to find a computationally efficient approximation of 
B a in the abstract setting of (40)-(41). 


Example 1 

The problem 0-0 fits in the abstract framework presented in this section 
when we assume that the state has regularity. We set W = L 2 (H), 

U = A = 1 — A, and T : H 2 {Bl) —> L 2 {dBl) is a trace operator, see 
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(44). Since A is a continuous isomorphism, assumptions |(Al)] and [(A2)| are 
both valid. The inner product on U a takes the form 


(u,v)u a = (Ku,v) + a {AM 1 Au,v) 

= uv ds + a (u — Au)(v — Av) dx 

Jon Jn 

= uvds + a D 2 u : D 2 v + 2X7u ■ Vu + uvdx, 

Jon Jn 

where D 2 u denotes the Hessian of u, and the last equality follows from the 
boundary condition du/d n = 0 imposed on H 2 (Q). The resulting precon¬ 
ditioner is the one that was used in the numerical experiments, detailed in 
Section [3| and it is spectrally equivalent to the preconditioner defined in 

dm 


Example 2 

Let U, W, and K be as in Example 1, but let us set A = —A. Now A has 
non-trivial kernel, consisting of the a.e. constant functions, and for constant 
u we have 

W Tu \\L 2 (dn) = J^ol~\\ u \\H 2 (n)- 

Since assumptions |(A1)| and |(A2)| are valid, the optimality system is still 
well-posed. In this case the inner product on U a is given by 


(u,v) 


u a = 


uvds + a / D 2 u:D 2 vdx. 


Idn 


Example 3 

Let us consider the “prototype” problem: 

mm{^||u-d||^ 2(f2) + |||/||i2 (n) } 

subject to 


—A u + u + f = 0 in 17, 

!"=() on dn. 

an 

Note that we here consider the case in which observation data is assumed 
to be available throughout the entire domain of the state equation. 

If the usual variational form of the PDE is used, i.e., 

0, w) H i(n) + (/, w) L 2 (U) =0, E H l (Q), (50) 
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then the control space equals L 2 (fl), whereas the dual space is The 

preconditioning strategy presented in this section is therefore not applicable. 

If instead we can assume i^ 2 (fi)-regularity, we can use the variational 
form 

{-Au + u,w) L 2 {Q ) + (f,w) L 2 in) =0, VmGl 2 (0). (51) 


Now, the control and dual spaces both equal L 2 ( Q). The methodology pre¬ 
sented in this section can thus be applied, and a robust preconditioner is 
obtained. Compared with the preconditioner for the problem with bound¬ 
ary observations only, see Section [5j equation (38), the only change is the 
replacement of Mg, in the (2, 2) block of B a with M. 

We remark that in m and [9], parameter-robust preconditioners were 
proposed for the “prototype” problem, using the standard variational formu¬ 
lation (50) of the PDE. Those methods do not require improved regularity 
for the state space. Instead, they require that observations are available 
throughout the computational domain. 


7 Eigenvalue analysis 


In Section 6] it was shown that the condition number of B a A a , with A a 
defined in (42b and B a defined in (49), can be bounded independently of a, 


as well as independently of the operators appearing in (40)-(41). Moreover, 
the numerical experiments indicate that the eigenvalues are contained in 
three intervals, independently of the regularization parameter a, see Figure 
[2] In this section we detail the structure of the spectrum of the precondi¬ 
tioned system considered in Section [6j and we obtain sharp estimates for 
the constants appearing in Theorem |6.1[ 


We consider self-adjoint linear operators A a and B a 


aM 

0 

M r 


aM 

0 

0 

0 

I< 

A 

and Hq, 1 = 

0 

K + aR 

0 

. M 

A 

0 


_ 0 

0 

a~ 1 M 


(52) 


where R is defined by 


r = Am- 1 a. 


(53) 


We assume that A : U —> W' and M : W -» W' are continuous operators, 
for some Hilbert spaces U and W. In addition we will make use of the 
following assumptions. 


(Bl) M is a self-adjoint and positive definite, 
(B2) I\ + R is positive definite, 

(B3) I\ is self-adjoint and positive semi-definite. 
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Assumptions | (B1) (B3)| ensure that B a is a self-adjoint and positive definite. 
In particular, assumptions (Bl) (B3) hold for A a as in (42), provided that 
the assumptions of Section [6] hold. For simplicity, we also assume that that 
A a and B a are finite-dimensional operators. 

Theorem 7.1 Let p, q, and r be the polynomials 

P( A) = 1 — A, q(\) = 1 +\p(\), r(A) =p- Xq(X). 

Let q\ < q 2 and r\ < r 2 < rj, be the roots of q and r, respectively. The 
spectrum of B a A a is contained within three intervals, determined by the 
roots of p and r, independently of a: 


sp(B a A a ) C [ 7 - 1 , 91 ] U [r-2,1] U [ 92 , 7 - 3 ]. 


(54) 


Consequently, the spectral condition number of B a A a is bounded, uniformly 
in a, 

Tn 

(55) 


7*3 


k(B a A a ) < — « 4.089. 
7-2 


If K has a nontrivial kernel, inequality (55) becomes an equality. 


Proof 

Consider the equivalent generalized eigenvalue problem 


otM 0 M' 
0 K A' 
MAO 



7' 



u 

= A 


w 



aM 0 0 

0 K + aR 0 
0 0 a~Hl 



7' 


u 


w 


(56) 


We show that (56) admits no nontrivial solutions unless A is as in 


Since M is a self-adjoint isomorphism, by assumption |(B1)[ we 
rewrite (56) as the three identities 


apf + w = 0, 
pKu + A'w — A aRu = 0, 

/ + M~ 1 Au — A a~ 1 w = 0. 

Assume that A is not contained within the three closed intervals of 


can 

(57) 

(58) 

(59) 


Then p A 0, and we can use (57) to eliminate / from 

0 = ap(f + M~ l Au — A a~ 1 iv) = apM~ l Au — (1 + A p)w 
= apM~ 1 Au — qw. 

Since q is nonzero, we can use ( |60| ) to eliminate w from 

0 = q{pKu + A’w — XaRu) = qpKu + a(p — \q)Ru 
= qpKu + rRu, 


(60) 


(61) 
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where the identity (53) was used. By assumption, pq and r are both nonzero. 
Moreover, it can be easily seen that pq and r have the same sign outside of 
the bounded intervals of (54). From assumptions (Bl) (B3)l we conclude 


that qpK + rR is a self-adjoint definite operator. Then (61) only admits 
trivial solutions, hence A can not be an eigenvalue of B a A a . 

The estimate (55) follows from (54), noting that |sp(i3 a „4 a )| C 
From (61) it can be seen that the roots of r are eigenvalues of B a A a if Ker K 
is nontrivial. 


Remark 

If A = (1 — A) : H 2 {£i) L 2 (Q)' , then R = A'M~ 1 A is characterized by a 

bilinear form &(■,•) as in ( |16[ ): 

(A , M~ 1 Au,v) = / AuAv + 2Vu • Vu + uv dx 

Jn 

= (u,v) H 2 (o)+ / Vu • Vu dx = b(u,v) 

Jn 

For discretizations Uh C U and Wh C W of A such that A(Uh) C M(Wh ), 
the discretization of b coincides with A h M^ 1 A^. This follows from an ar¬ 
gument similar to that in the proof of Lemma |4.6[ , and as a consequence, 
Theorem |7.1| can be applied to the preconditioned discrete systems consid¬ 
ered in Section [3j 

8 Discussion 

Previously, parameter robust preconditioners for PDE-constrained optimiza¬ 
tion problems have been successfully developed, provided that observation 
data is available throughout the entire domain of the state equation. For 
many important inverse problems, arising in industry and science, this is 
an unrealistic requirement. On the contrary, observation data will typically 
only be available in subregions, of the domain of the state variable, or at 
the boundary of this domain. We have therefore explored the possibility for 
also constructing robust preconditioners for PDE-constrained optimization 
problems with limited observation data. 

For an elliptic control problem, with boundary observations only, we have 
developed a regularization robust preconditioner for the associated KKT sys¬ 
tem. Consequently, the number of Minres iterations required to solve the 
problem is bounded independently of both regularization parameter a and 
the mesh size h. In order to achieve this, it was necessary to write the ellip¬ 
tic state equation on a non-standard, and non-self-adjoint, variational form. 
If this approach is employed, then the control and the Lagrange multiplier 
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will belong to the same Hilbert space, which leads to extra regularity re¬ 
quirements for the state. This fact makes it possible to construct parameter 
weighted metrics such that the constants appearing in the Brezzi conditions, 
as well as the constants in the inequalities expressing the boundedness of the 
KKT system, are independent of a and h. Consequently, the spectrum of 
the preconditioned KKT system is uniformly bounded with respect to a and 
h, which is ideal for the Minres scheme. These properties were illuminated 
through a series of numerical experiments, and the preconditioned Minres 
scheme handled our model problem excellently. 

The use of a non-self-adjoint form of the elliptic state equation leads 
to additional challenges for constructing discretization schemes and suitable 
multigrid methods. More specifically, it becomes necessary to implement a 
FE space approximating H 2 . We accomplished this by a C 1 discretization 
that is conforming in H 2 . The method employed does, however, have strong 
restrictions on the mesh, which seemingly must be composed of rectangles. 


References 

[1] F. Brezzi. On the existence, uniqueness and approximation of saddle- 
point problems arising from Lagrangian multipliers. RAIRO numerical 
analysis , 8:129-151, 1974. 

[2] P. Grisvard. Elliptic problems in nonsmooth domains. Pitman, Boston, 
1985. 

[3] A. Gunnel, R. Herzog, and E. Sachs. A note on preconditioners and 
scalar products in Krylov subspace methods for self-adjoint problems in 
Hilbert space. Electronic Transactions on Numerical Analysis , 41:13- 
20, 2014. 

[4] W. Hackbusch. Elliptic differential equations. Theory and numerical 
treatment. Springer-Verlag, 1992. 

[5] K. A. Mardal and R. Winther. Preconditioning discretizations of sys¬ 
tems of partial differential equations. Numerical Linear Algebra with 
Applications , 18(l):l-40, 2011. 

[6] B. F. Nielsen and K. A. Mardal. Efficient preconditioners for optimality 
systems arising in connection with inverse problems. SIAM Journal on 
Control and Optimization , 48(8), 2010. 

[7] B. F. Nielsen and K. A. Mardal. Analysis of the Minimal Residual 
Method applied to ill-posed optimality systems. SIAM Journal on Sci¬ 
entific Computing , 35(2):A785-A814, 2013. 


26 



[8] J. W. Pearson, M. Stoll, and A. J. Wathen. Regularization-robust pre¬ 
conditioners for time-dependent PDE-constrained optimization prob¬ 
lems. SIAM Journal on Matrix Analysis and Applications , 33:1126- 
1152, 2012. 

[9] J. W. Pearson and A. J. Wathen. A new approximation of the Schur 
complement in preconditioners for PDE-constrained optimization. Nu¬ 
merical Linear Algebra with Applications, 19:816-829, 2012. 

[10] J. Schoberl and W. Zulehner. Symmetric indefinite preconditioners for 
saddle point problems with applications to PDE-constrained optimiza¬ 
tion problems. SIAM Journal on Matrix Analysis and Applications, 
29(3):752-773, 2007. 

[11] X. Zhang. Multilevel schwarz methods for the biharmonic dirichlet 
problem. SIAM Journal on Scientific Computing, 15(3):621-644, 1994. 

[12] W. Zulehner. Nonstandard norms and robust estimates for saddle point 
problems. SIAM Journal on Matrix Analysis and Applications, 32:536- 
560, 2011. 


27 



